import numpy as np
from scipy import signal,misc,ndimage
from skimage import filters,feature,img_as_float
from skimage.io import imread
from skimage.color import rgb2gray
from PIL import Image,ImageFilter
import matplotlib.pylab as plt
from skimage.filters import laplace
from scipy.signal import convolve2d
from scipy.ndimage import gaussian_filter
from PIL.ImageFilter import (FIND_EDGES,EDGE_ENHANCE,EDGE_ENHANCE_MORE)
from skimage.transform import pyramid_gaussian
from skimage.transform import pyramid_reduce
from skimage.transform import pyramid_expand,resize
from skimage.util import random_noise





plt.rcParams['font.sans-serif']=['simhei']
plt.rcParams['font.family']='sans-serif'


# 5.1.1
# 使用卷积核来计算梯度，以灰度图输入，绘制图像像素值和梯度分量的x分量随y的变化而变化
def plot_image(image,title):
    plt.imshow(image)
    plt.title(title,size=20)
    plt.axis('off')

ker_x=[[-1,1]]
ker_y=[[-1],[1]]
im=rgb2gray(imread('zebra.png'))
im_x=signal.convolve2d(im,ker_x,mode='same')
im_y=signal.convolve2d(im,ker_y,mode='same')
im_mag=np.sqrt(im_x**2+im_y**2)
im_dir=np.arctan(im_y/im_x)
plt.gray()
plt.figure(figsize=(30,20))
plt.subplot(231),plot_image(im,'原始图像')
plt.subplot(232),plot_image(im_x,'grad_x')
plt.subplot(233),plot_image(im_y,'grad_y')
plt.subplot(234),plot_image(im_mag,'||grad||')
plt.subplot(235),plot_image(im_dir,r'$\theta$')
plt.subplot(236)
plt.plot(range(im.shape[1]),im[0,:],'b-',label=r'$f(x,y)|_{x=0}$',linewidth=5)
plt.plot(range(im.shape[1]),im_x[0,:],'r-',label=r'$grad_x(f(x,y))|_{x=0}$')
plt.title(r'$grad_x (f(x,y))|_{x=0}$',size=30)
plt.legend(prop={'size':20})
plt.show()


# 5.1.2
# 使用核卷积来计算图像的拉普拉斯算子
ker_laplacian=[[0,-4,0],[-1,4,-1],[0,-1,0]]
im=rgb2gray(imread('zebra.png'))
im1=np.clip(signal.convolve2d(im,ker_laplacian,mode='same'),0,1)
plt.gray()
plt.figure(figsize=(20,10))
plt.subplot(121),plot_image(im,'原始图像')
plt.subplot(122),plot_image(im1,'拉普拉斯算子的卷积')
plt.show()


#5.1.3
sigma=1
im=im+random_noise(im,var=sigma**2)
plt.imshow(im)
plt.show()


# 5.2.1
# 应用拉普拉斯滤波器，得到的结果与原始图像叠加，使用scikit-image filters模块的laplace（）函数实现
im=rgb2gray(imread('panda.png'))
im1=np.clip(laplace(im)+im,0,1)
plt.figure(figsize=(20,30))
plt.subplot(211),plot_image(im,'原始图像')
plt.subplot(212),plot_image(im1,'锐化图像')
plt.tight_layout()
plt.show()


# 5.2.2
# 使用SciPy的ndimage模块实现范蕊花掩膜
def rgb2gray(im):
    return np.clip(0.2989*im[...,0]+0.5870*im[...,1]+0.1140*im[...,2],0,1)
im=rgb2gray(img_as_float(imread('man.png')))
im_blurred=ndimage.gaussian_filter(im,5)
im_detail=np.clip(im-im_blurred,0,1)
plt.gray()
fig, axes=plt.subplots(nrows=2,ncols=3,sharex=True,sharey=True,figsize=(15,15))
axes=axes.ravel()
axes[0].set_title('原始图像',size=15),axes[0].imshow(im)
axes[1].set_title('模糊图像,sigma=5',size=15)
axes[1].imshow(im_blurred)
axes[2].set_title('细节图像',size=15),axes[2].imshow(im_detail)
alpha=[1,5,10]
for i in range(3):
    im_sharp=np.clip(im+alpha[i]*im_detail,0,1)
    axes[3+i].imshow(im_sharp),axes[3+i].set_title('锐化图像,alpha='+str(alpha[i]),size=15)
for ax in axes:
    ax.axis('off')
fig.tight_layout()
plt.show()


# 5.3.2
# 利用函数sobel-h（）、sobel-y（）、sobel分别查找水平/垂直边缘，并使用Sobel算计计算梯度大小
im=rgb2gray(imread('tajmahal.jpg'))
plt.gray()
plt.figure(figsize=(20,18))
plt.subplot(2,2,1)
plot_image(im,'原始图像')
plt.subplot(2,2,2)
edges_x=filters.sobel_h(im)
plot_image(edges_x,'sobel_x')
plt.subplot(2,2,3)
edges_y=filters.sobel_v(im)
plot_image(edges_y,'sobel_y')
plt.subplot(2,2,4)
edges=filters.sobel(im)
plot_image(edges,'sobel')
plt.subplots_adjust(wspace=0.1,hspace=0.1)
plt.show()


# 5.3.3
# 对图像应用不同的边缘检测滤波器，得到梯度大小
im=rgb2gray(imread('vic.png'))
plt.gray()
plt.figure(figsize=(20,24))
plt.subplot(3,2,1)
plot_image(im,'原始图像')
edges=filters.roberts(im)
plt.subplot(3,2,2)
plot_image(edges,'应用Roberts滤波器的图像')
edges=filters.scharr(im)
plt.subplot(3,2,3)
plot_image(edges,'应用Roberts滤波器的图像')
edges=filters.sobel(im)
plt.subplot(3,2,4)
plot_image(edges,'应用Roberts滤波器的图像')
edges=filters.prewitt(im)
plt.subplot(3,2,5)
plot_image(edges,'应用Roberts滤波器的图像')
edges=np.clip(filters.laplace(im),0,1)
plt.subplot(3,2,6)
plot_image(edges,'应用Roberts滤波器的图像')
plt.subplots_adjust(wspace=0.1,hspace=0.1)
plt.show()


# 5.3.4
# Canny边缘检测器的实现方式
im=rgb2gray(imread('tiger.jpg'))
im=ndimage.gaussian_filter(im,4)
im+=0.05*np.random.random(im.shape)
edges1=feature.canny(im)
edges2=feature.canny(im,sigma=3)
fig,(axes1,axes2,axes3)=plt.subplots(nrows=1,ncols=3,figsize=(30,12),sharex=True,sharey=True)
axes1.imshow(im,cmap=plt.cm.gray),axes1.axis('off')
axes1.set_title('噪声图像',fontsize=50)
axes2.imshow(edges1,cmap=plt.cm.gray),axes2.axis('off')
axes2.set_title('Canny滤波器,$\sigma=1$',fontsize=50)
axes3.imshow(edges2,cmap=plt.cm.gray),axes3.axis('off')
axes3.set_title('Canny滤波器,$\sigma=3$',fontsize=50)
fig.tight_layout()
plt.show()


# 5.3.5
# 计算LoG滤波器核相应最佳DoG近似（已知值），将其应用在相同的输入图像上
def plot_kernel(kernel,s,name):
    plt.imshow(kernel,cmap='YlOrRd')
    ax=plt.gca()
    ax.set_xticks(np.arange(-0.5,kernel.shape[0],2.5))
    ax.set_yticks(np.arange(-0.5,kernel.shape[1],2.5))
    plt.colorbar()


def LOG(k=12,s=3):
    n=2*k+1
    kernel=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            kernel[i,j]=-(1-((i-k)**2+(j-k)**2)/(2.*s**2))*np.exp(-((i-k)**2+(j-k)**2)/(2.*s**2))/(np.pi*s**4)
    kernel=np.round(kernel/np.sqrt((kernel**2).sum()),3)
    return kernel

def DOG(k=12,s=3):
    n=2*k+1
    s1,s2=s*np.sqrt(2),s/np.sqrt(2)
    kernel=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            kernel[i,j]=np.exp(-((i-k)**2+(j-k)**2)/(2.*s1**2))/(2*np.pi*s1**2)-np.exp(-((i-k)**2+(j-k)**2)/(2.*s2**2))/(2*np.pi*s2**2)
    kernel=np.round(kernel/np.sqrt((kernel**2).sum()),3)
    return kernel
s=3
img=rgb2gray(imread('face.png'))
kernel=LOG()
outimg=convolve2d(img,kernel)
plt.figure(figsize=(20,20))
plt.subplot(221)
plt.title('LOG kernel',size=20),plot_kernel(kernel,s,'DOG')
plt.subplot(222)
plt.title('output image with LOG',size=20)
plt.imshow(np.clip(outimg,0,1),cmap='gray')
kernel=DOG()
outimg=convolve2d(img,DOG())
plt.subplot(223)
plt.title('DOG kernel',size=20),plot_kernel(kernel,s,'DOG')
plt.subplot(224)
plt.title('output image with DOG',size=20)
plt.imshow(np.clip(outimg,0,1),cmap='gray')
plt.show()

# 基于SciPy ndimage 模块的LoG滤波器
img=rgb2gray(imread('zebra.png'))
fig=plt.figure(figsize=(25,15))
plt.gray()
for sigma in range(1,10):
    plt.subplot(3,3,sigma)
    img_log=ndimage.gaussian_laplace(img,sigma=sigma)
    plt.imshow(np.clip(img_log,0,1))
    plt.axis('off')
    plt.title('LoG with sigma='+str(sigma),size=20)
plt.show()


# 5.3.6
# 利用零点交叉检测图像的边缘
fig=plt.figure(figsize=(25,15))
plt.gray()
for sigma in range(2,10,2):
    plt.subplot(2,2,sigma/2)
    result=ndimage.gaussian_laplace(img,sigma=sigma)
    plt.imshow(result)
    plt.axis('off')
    plt.title('LoG with zero-crossing,sigma='+str(sigma),size=20)
plt.show()


# 5.3.7
# 基于PIL发现和增强边缘
im=Image.open('wordtrade.jpg')
plt.figure(figsize=(18,25))
plt.subplot(2,2,1)
plot_image(im,'original (UMBC library)')
i=2
for f in (FIND_EDGES,EDGE_ENHANCE,EDGE_ENHANCE_MORE):
    plt.subplot(2,2,i)
    im1=im.filter(f)
    plot_image(im,str(f))
    i+=1
plt.show()


# 5.4.1
# scikit-image transform pyramid 模块的高斯金字塔
image= imread('lena.jpg')
nrows,ncols=image.shape[:2]
pyramid=tuple(pyramid_gaussian(image,downscale=2))
plt.figure(figsize=(20,5))
i,n=1,len(pyramid)
for p in pyramid:
    plt.subplot(1,n,i),plt.imshow(p)
    plt.title(str(p.shape[0])+'x'+str(p.shape[1])),plt.axis('off')
    i+=1
plt.subplot('Gaussian Pyramid',size=30)
plt.show()
compos_image=np.zeros((nrows,ncols+ncols//2,3),dtype=np.double)
compos_image[:nrows,:ncols,:]=pyramid[0]
i_row=0
for p in pyramid[1:]:
    nrows,ncols=p.shape[:2]
    compos_image[i_row:i_row+nrows,cols:cols+ncols]=p
    i_row+=nrows
fig,axes=plt.subplots(figsize=(20,20))
axes.imshow(compos_image)
plt.show()


# 5.4.3
# 构造高斯金字塔
def get_gaussian_pyramid(image):
    gaussian_pyramid=[]
    return gaussian_pyramid

def get_laplacian_pyramid(gaussian_pyramid):
    laplacian_pyramid=[]
    return laplacian_pyramid

image=imread('bunny.png')
gaussian_pyramid=get_gaussian_pyramid(image)
laplacian_pyramid=get_laplacian_pyramid(gaussian_pyramid)
w,h=20,12
p=gaussian_pyramid[0]
for i in range(3):
    plt.figure(figsize=(w,h))
    p=gaussian_pyramid[i]
    plt.imshow(p)
plt.title(str(p.shape[0])+'x'+str(p.shape[1]),size=20)
plt.axis('off')
w,h=w/2,h/2
plt.show()

w,h=10,6
for i in range(1,4):
    plt.figure(figsize=(w,h))
    p=laplacian_pyramid[i]
    plt.imshow(rgb2gray(p),cmap='gray')
plt.title(str(p.shape[0])+'x'+str(p.shape[1]),size=20)
plt.zxis('off')
w,h=w/2,h/2
plt.show()


# 5.4.4
# 仅通过拉普拉斯金字塔重建图像
def reconstruct_image_from_laplacian_pyramid(pyramid):
    i=len(pyramid)-2
    prev=pyramid[i+1]
    plt.figure(figsize=(20,20))
    j=1
    while i>=0:
        prev=resize(pyramid_expand(prev,upscale=2),pyramid[i].shape)
        im=np.clip(pyramid[i]+prev,0,1)
        plt.subplot(3,3,j)
        plt.imshow(im)
        plt.title('Level='+str(j)+''+str(im.shape[0])+'x'+str(im.shape[1]),size=20)
        prev=im
        i-=1
        j+=1
        plt.subplot(3,3,j)
        plt.imshow(image)
        plt.title('Original image'+''+str(image.shape[0])+'x'+str(image.shape[1]),size=20)
        plt.show()
        return im

image=img_as_float(imread('deer.png')[...,:3])
pyramid=get_laplacian_pyramid(get_gaussian_pyramid(image))
im=reconstruct_image_from_laplacian_pyramid(pyramid)


# 5.4.5
# 基于金字塔的图像融合
A=img_as_float(imread('ronaldo.png')[...,:3])
B=img_as_float(imread('messi.png')[...,:3])
M=img_as_float(imread('deer.png')[...,:3])

pyramidA=get_laplacian_pyramid(get_gaussian_pyramid(A))
pyramidB=get_laplacian_pyramid(get_gaussian_pyramid(B))
pyramidM=get_gaussian_pyramid(M)

pyramidC=[]
for i in range(len(pyramidM)):
    im=pyramidM[i]*pyramidA[i]+(1-pyramidM[i])*pyramidB[i]
    pyramidC.append(im)

I=reconstruct_image_from_laplacian_pyramid(pyramidC)